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ABSTRACT 

Two crystal growth processes being considered for Spacelab 
experiments were studied in this program in order to permit antici- 
pation and understanding of phenomena not ordinarily encountered 
on earth. Computer calculations were performed on transport pro- 
cesses in floating zone melting and on growth of a crystal from 
solution in a spacecraft environment. Experiments intended to 
simulate solution growth at micro accelerations were performed. 

Surface-tension driven convection was found to be vigorous 
in floating zone melting of silicon. It is relatively unaffected by 
zone movement and by buoyancy effects even at earth's gravity. With 
very effective heat shielding, a laminar cell resembling a smoke ring 
is set up in the top half of the zone, and another in the lower half. 
Mass transfer is slow between these cells. It is also less effective 
along the center of the zone than at the periphery, leading to a 
significant radial variation in composition in the resulting crystal, 
unless very low growth rates are employed. Less effective heat 
shielding would produce secondary cells in the melt. These would 
probably become oscillatory and cause compositional striations parallel 


to the growth interface. Without heat shielding the convection would 
almost certainly be turbulent. Turbulence would result in a fairly 
uniform composition radially, but would cause very fine compositional 
striations, probably on a scale too small to be detected by ordinary 
means. Electrical and optical properties would probably be degraded, 
however. In good thermal conductors, the surface driven convection 
has only a small effect on heat transfer. In poor thermal conductors, 
the influence is large. As a consequence, one might expect a twitch- 
ing or oscillatory convection at the melt surface when a poor con- 
ductor is heated by radiation along the entire melt. 

Both our transient numerical calculations and previous steady 
state theory predict that a vortex-like ring should be established 
about a suspended sphere slowly growing from a solution at low g. 
However our low g simulation experiments produce a stream which 
separates from the sphere and meanders about the solution, only slowly 
mixing in. The solution soon became very inhomogeneous. The con- 
vection was insensitive to mechanical perturbations. These results 
are to be contrasted with the usual solution-growth system in which 
the solution density depends strongly on composition. The stream 
which separates from the crystal in this case rapidly mixes with the 
bulk solution. Small mechanical perturbations produce "explosions" 
in this convection— a sudden increase in convection, with turbulent 
eddies, which then decrease again to steady state. Perhaps the most 
encouraging result was that inclusion-free faceted cyrstals were pro- 
duced under simulated low g conditions. 


iii 



Two papers on this work were delivered at the May 17 San Diego 
meeting of the Western Section of the American Association of Crystal 
Growth. A talk on the floating zone studies was given at the June 
1974 meeting of the American Institute of Chemical Engineers in 
Pittsburgh. One paper (J. Crystal Growth 21^ (1974) 182) was published 
and three others are in preparation. 
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SUMMARY 


Surface-tension driven convection in floating zone melting 
was studied by computer solution of uhe partial differential equa- 
tions for steady-state momentum, heat and mass transfer. Finite 
difference methods were employed. Two situations were investigated — 
a parabolic temperature profile along the melt surface, as might 
result from radiant heating; and a line heat source as with electron 
beam h'.ating. Although tha calculations were performed using 
!imensionless parameters, the conclusions below are expressed in 
terms of a 1 cm diameter x 1 cm long molten zone of silicon, except 
where noted. 

1. The flow patterns for both methods of heating are similar, 
but with radiant heating the highest fluid velocities are at the melt 
periphery near the crystal interface, and with electron beam heating 
producing the highest velocities near the middle of the zone. 

2. For a temperature difference of 0.005°C between the inter- 
face and the middle of the zone at the periphery, a single roll cell 
is formed in the top half of the zone, with its mirror image in the 
bottom half. The maximum melt velocity is about 3 cm/min. Such a 
small temperature difference could only be produced by use of very 
effective heat shields or after-heaters. 

3. At about 0.1°C temperature difference, secondary cells are 
produced. Although our steady state calculations could not show it, 
we suspect that the flow would become oscillatory and eventually 
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turbulent as this temperature difference is increased (by reducing the 
heat shielding) . 

4. For good thermal conductors such as silicon, the fluid 
motion has little influence on heat transfer and interface position. 
The reverse is true for poor conductors. 

5. Normal zone travel rates have little influence on the 
surface-driven flow. For example, with a temperature difference of 
0.005°C, a zone travel rate of 10 cm/hour noticeably perturbs the 
streamlines only in the vicinity of the solid-liquid in r arfaces. 

6. At small temperature gradients, the primary resistance to 
mass transfer is in the middle of the zone, unlike boundary layer 
flow in which the mass transfer resistance is concentrated near the 
solid-liquid interfaces. 

7. Significant radial inhomogeneities in the resulting 
crystal would be expected when the convection is relatively slow. For 
example, a 10% radial variation would result from a zone travel rate 
of 1 cm/hour, a distribution coefficient of 0.1, and a temperature 
difference of 0.05°C along the melt surface. At 4 cm/hour zone 
travel rate there would be a 40% impurity concentration gradient. 

8. Buoyancy-driven flow at earth's gravity is negligible com- 
pared to surface driven flow for molten silicon. 

9. The vigor of surface driven flow varies considerably from 
material to material, increasing as the dependence of surface tension 
on temperature increases, as density increases, and as the viscosity 
of the melt decreases. Since temperature gradients generally increase 



as thermal conductivity decreases, one would also expect poor thermal 
conductors to have greater surface* driven flow, other properties being 
equal. 

Theoretical predictions for solution growth at low accelerations 

are: 

1. The flow will not be of the boundary-layer type. The 
perturbed region will extend into the solution for distances on the 
order of the crystal size. 

2. Long times will be required for convection to influence the 
mass transfer and for steady state to be reached. 

Several solute-solvent systems have been found which have a 
density only weakly dependent on composition. These are useful for 
simulation of low acceleration solution growth and dissolution because 
they yield Grashof (Gr) numbers on earth of the same magnitude as for 
normal crystal growth systems in a space laboratory. Optical and 
photographic techniques have been developed for studying growth and 
dissolution. Preliminary experiments have revealed that: 

1. High quality crystals can be grown from solution under low 
Grashof (Gr) number conditions, as exemplified by thymol growing from 
ethyl carbonate solutions. 

2. At low Gr the stream formed during dissolution mixes much 
more slowly with the bulk solution than at high Gr. 

3. The convective flow pattern at low Gr is much less sensitive 
to perturbations than at high Gr. 


xiv 




The calculational portions of this program have ended, pending 
additional funding. The experiments will continue with present funds 
until December. 




1 


} 


* 



I. INTRODUCTION 


This program consisted of three parts: 

1. Computer calculations of surface tension driven flow in 

floating zone melting at zero g. 

2. Computer calculations of transient free convection about 
a suspended sphere at very low g. 

3. Experimental simulation of solution crystal growth in a 
space environment. 

Part 1 was initiated in late June 1973, Part 2 in August 1973 and Part 
3 in December 1973. Without additional funding. Parts 1 and 2 have 
ended. 

One of the most promising space processing applications is 
floating zone melting. On earth the crystal diameter that can be 
produced is limited by the size of the zone that can be supported by 
the surface tension of the melt in opposition to gravity. In space 
any size could be produced. For many materials the space vacuum could 
advantageously be utilized as an environment for the float zoning. 
While the absence of earth's gravity assures that buoyancy-driven free 
convection will be negligible in space, one must not assume that there 
will be no convection in the zone. Since the temperature will vary 
along the gas-melt surface of the zone, the surface tension will vary 
along the surface of the zone, normally increasing as one moves from 


! 
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the center toward the solid. This surface tension gradient will 
develop a flow field, as the system tries to minimize its surface 
energy via movement of the melt surface from low surface tension 
regions to high surface tension regions. The question is, how 
vigorous is this flow and what influence does it have on the crystal 
growth? One cannot do experiments on earth because of interference 
with buoyancy effects. Even in space only limited experiments could 
be done on opaque materials such as silicon. However, the governing 
differential equations can be solved on earth to provide a fairly 
complete picture of surface-driven flow phenomenon. We have solved 
some oi the interesting problems here. 

Before beginning this program there was some concern on our 
part about the advisability of growing crystals from solution in 
space, because of the grc?t tendency to occlude solution during growth. 
However, our experiments have shown that inclusion-free faceted crystals 
;an be produced in a low-convection environment. Since there will be 
some acceleration in Spacelab, free convection will not be completely 
absent. The convection will probably not be steady state, however, 
because of the long transients predicted and the variation in 
a -'-deration in an orbiting laboratory. Our objective is to determine 
the convection patterns to be expected and their influence on crystal 
growth phenomena. Both simulation experiments and transport calcula- 
tions have been performed. Solutions with a density nearly independent 
of composition are required for the experiments. 



Each of the following three sections of this report covers 
one of the parts described abov-j. These sections are independent 
of one another so that it is not necessary to read one section in 
order to understand another. Possible future research is discussed 


in the final section 
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II. SURFACE TENSION DRIVEN FLOW IN FLOATING 
ZONE MELTING AT ZERO G 


(Chong E. Chang) 


Definition of Symbols 


a 

b 


c 


P 

C 


C 


ct 


i 
i#t 


»b 


D 


g 


g 


e 


Gr h 


Gr 


m 


i 

j 

k 


Radius of the molten zone, cm 

Distance of the top surface of the solid from the 
center of the zone, cm 

Specific heat, cal/g°K 

Dimensionless impurity concentration in the liquid, 
defined by m/m^ 

Dimensionless impurity concentration in the feed 
crystal, defined by m cfc /m t 

Dimensionless impurity concentration in the liquid 

at the crystallizing interface, defined by m. ,/m 

1 ,b t 

Dimensionless impurity concentration in the liquid 

at the melting interface, defined by m. 

i,t t 

2 

Binary molecular diffusion coefficient, cm /sec 

2 

Acceleration, cm/sec 

2 

Earth's gravitational field, 980 cm/sec 

Grashof number for heat transfer (~ ratio of the 

temperature difference buoyancy force to the 

viscous force) defined by g$(T -T 

o m 

Grashof number for mass transfer (~ ratio of the 
composition buoyancy force to the viscous force) 
defined by gam^a-Vv 2 

Integer denoting grid station in axial direction, 
i - 1,2,3,... 

Integer denoting grid station in radial direction, 
j m 1,2,3,... 

Thermal conductivity, cal/ cm sec°C 
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s 


k Inter facial distribution coefficient (ratio of the 

w impurity concentration in the crystal to that in 
the liquid at the interface) 

i One half of the liquid zone length, cm 

3 

m Impurity concentration in liquid, 7/cm 


m. 


m 


cb 


m 


ct 


m 


i,b 


m 


i*t 


M' 


Nu 


P 

P o 

P 

Pr 


r 


Inpurity concentration the feed solid would have 


if it were melted, g/cm 3 

Impurity concentration in 

Impurity concentration in 

Impurity concentration in 
lizing interface, g/cm^ 

Impurity concentration in 
interface, g/cm 3 


3 

the crystal, g/cm 

3 

the feed crystal, g/cm 
the liquid at the crystal- 

the liquid at the melting 


Dimensionless surface tension parameter at the free 

liquid surface (~ ratio of surface tension force to 

viscous force) defined by pa (T -T ) a/y^ 

o m 

Nusselt number (dimensionless temperature gradient, 
or ratio of total heat transport to heat conduction) 
defined by | 90 / 3 z] + ^ a 

2 • 

Dynamic pressure in fluid dyne/cm , P-p Q 

2 

Static pressure in fluid, dyne/cm , defined by 
p f gU-z) 

2 

Total pressure in fluid, dyne/cm 

Prandtl number, yC /k (ratio of the momentum dif- 
fusivity to the thlrmal diffusivity) 

Radial coordinate, measured from the center of the 
zone, cm 


R Dimensionless radial coordinate, measured from the 
center of the zone, r/a 

Ar Dimensionless radial distance between grid points; 


Sc Schmidt number (ratio of the momentum diffusivity 
to the molecular diffusion coefficient) , y/pD 
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Sh Sherwood number (dimensionless concentration 
gradient), defined by |3c/3z| +J ^, a 

o 

T Temperature, K 

T Temperature of the middle of liquid zone at the 

periphery, °K 

T c Temperature of the surroundings, °K 

o 

T Interfacial temperature (melting point) , K 

m 

AT Either (T Q -T m ) if 0 is defined by the Equation (A- 
19) or (T m -T c ) if 0 is defined by the Equation (A- 
31), °K 

v c Zone travel rate, cm/sec 

v^ Radial velocity component of the fluid, cm/sec 

v Axial velocity component of the fluid, cm/sec 

z 

V Dimensionless zone travel rate, v a/v 

c c 

V Dimensionless radial velocity component of the 

r fluid, v r a/V 

V Dimensionless axial velocity component of the 

Z fluid, v a/V 

Z 

V . Dimensionless interfacial flow rate, V (p /p,) 
ci c c c 

z Axial coordinate, measured from the center of the 
zone , cm 

Z Dimensionless axial coordinate, measured from the 
center of the zone, z/a 

Az Dimensionless axial distance between grid points 

a (3p/3m) T( p/P; volume expansion coefficient due to a 
concentration change in fluid, l/(g/cm^) 

0 (3p/3T) m ,p/P; volume expansion coefficient due to a 

temperature change in fluid, l/°c 

Y Surface tension between melt and vapor, dyne/cm 

e # Emissivity of thermal radiation 
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y 

V 

p 


' f 
0 


♦ 


U) 


Dimensionless temperature, (T-T )/At 

m 

Dimensionless temperature of the surroundings, 

(T -T )/AT 
c m 

Viscosity of the fluid, g/cm sec 

2 

Kinematic viscosity, \i/ p, cm /sec 
Density of the fluid, g/cm^ 

3 

Density of the crystal, g/cm 

3 

Density of the fluid at melting point, g/cm 

-5 2 

Stefan-Boltzmann constant, 5.668 x 10 erg/cm 
0 4 
sec K 

Temperature coefficient of the surface tension, 
3y/3T, dynes/cm°K 

Dimensionless stream function, defined by 
Equations (A-25) and (A-26) in Appendix A. 

Dimensionless vorticity, defined by Equations (A-21) 
and (A-27) in Appendix A. 


A. Introduction 

Floating zone melting (1) is a common method of growing single 
crystals, especially for high melting materials because the melt is 
not in contact with a crucible. Although free convection due to 
density variations will not occur at zero g, temperature variations 
along the free melt surface will cause surface tension variations 
which generate convection. Surface-driven flow is often called the 
"Marangoni effect". 

We consider a liquid zone suspended between two cylindrical 
solids by its own surface tension, as in F gure 1. Both the solid and 
the liquid are surrounded by a gas or vacuum. The liquid zone is 


t 
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generated by a power source , typically a radiant heater, an electron 
bean, or by induction from a high frequency coil. (Induction heating 
is generally considered too inefficient and bulky to be employed for 
space processing, although it is used extensively on earth for con- 
ductors.) The zone is made to move by slowly moving the heater, 


SOLID 


LIQUID 


Z = l 


GAS OR 
VACUUM 


Z =-Z 


r * a 


SOLID 


Figure 1. Floating Zone Melting, Geometry, Coordinate 
System and Fluid Velocity Components 


\ ' 


resulting in continuous solidification at one liquid/solid interface 
and melting at the other. In order to determine the influence of 
surface driven flow on zero melting, the partial differential 
equations for momentum and heat transfer were considered. These 
are shown in the appendix. Solutions were obtained num> rlcally usir j 
the computer for two situations — a parabolic temperature profile 
along the melt-vapor surface and for a ring heat source at the center 
of the zone. The parabolic profile corresponds approximately to 
radiant heating, and the ring source corresponds to electron beam 
heating. The mass transfer equation coupled with the computed flow 
fields were employed to find the impurity concentration fields in the 
melt and in the growing crystal at steady state. Most of the calcu- 
lations were done for a 1 cm diameter silicon rod, although the 
results are expressed in dimensionless form so that they also represent 
other possibilities. 

B. Model 

The coordinate system and geometry are shown schematically in 
Figure 1. The change in modes of heating can be easily taken into 
account by simply changing the thermal boundary conditions. Since 
we consider only steady state, the Bulerian specification of the flow 
field seems to be convenient. That is, the interfaces and heater are 
taken to be fixed so that material continuously moves through the zone 
and the heater. This is equivalent to suction at the freezing inter- 
face and blowing at the melting interface (4) , 


Assumptions used in the analysis are as follows t 

a) Axisymmetrical flow) i.e., no angular components. (Violated 
if crystal rotated or if flow is not steady.) 

b) Incompressible liquid. (Variation of p negligible.) 

c) Constant properties of the liquid zone, i.e., p, y, k, D, 

Cp, e, all constant. 

d) Flanar liquid/solid interfaces. 

e) ohape of the liquid zone is circular cylinder. 

f) Symmetrical heating of the zone. 

g) Negligible viscous dissipation (low Mach n .: .*). 

h) Steady state. 

i) Growth and melting interface kinetics are sufficiently rapid 
that the liquid/solid interfaces are at the melting point of 
the material. 

j) Constant distribution coefficient k . 

o 

The mass conservation, Navier-Stokes , mass transfer equations, 
and boundary conditions are given in Appendix A, along with a descrip- 
tion of the methods of solution. Silicon was taken as a model material 
for numerical calculation of flow, temperature, and impurity (dopant) 
concentration fields. For silicon, the melting point is 1410°C, the 

Prandtl number is 0.023, the Schmidt number is about 5, M'/a(T -T ) - 

o m 

2 

p 0/\i is 14,000, and the emiBsivity is assumed to be 0.3 for both the 
melt ar.d the solid. 


Results of the calculations are summarized in Table I 
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C. Streamlines for a Parabolic Temperature Profile at 
the Free Melt Surface (Radiation Heating) _____ 

For simplicity, a parabolic temperature profile along the melc. 
surface was assumed for our initial calculations. This enabled us to 
solve the momentum equations without simultaneously solving the heat 
transfer equations. A parabolic profile, with a maximum at the center, 
seems reasonable for radiant heating. The dimensionless surface 
tension parameter M' increases as the radius of the zone a and the 
temperature variation along the melt surface increases. In Figures 
2 through 5, the streamlines for parabolic temperature profiles are 
shown to illustrate the effect of increased values of M (35, 350, 

3500 and 7000) . Donut-shaped vortex cells were formed. With M' = 

35 and 350 only two cells are generated, and the centers of the 
vortices move closer toward the liquid/solid interfaces as M* in- 
creases, as shown in Figures 2 and 3. With M' = 3500, however, which 

corresponds to a condition of a = 0.5 cm and (T -T ) = 0.5°C for 

o m 

silicon, secondary vortex cells were induced behind the primary 
vortex, as shown in Figure 4. As the value of M' was further in- 
creased up to 7000, third and fourth vortices were produced, as shown 
in Figure 5. These multiple vortices are probably indicative of 
oscillations and incipient turbulence, which cannot be calculated in 
a steady state analysis. Oscillations are frequently found with free 
convection in enclosed cavities between laminar and fully turbulent 

flow. The maximum velocities for 1 cm diameter silicon with (T -T ) » 

o m 

1.0°C are calculated to be 2 cm/sec. 
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Figure 2. Computed dimensionless streamlines ^ for surface 

driven flow in a floating zone at zero gravity with 
a parabolic temperature profile on the free liquid 
surface with M' ■ 35 and v Q » 0. This would be 
expected, for example, for silicon with (T 0 -T ro ) - 
0.005C and a - 0.5 cm. 
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Figure 3. Computed dimensionless streamlines ^ for surface tension 
driven flow in a floating zone at zero gravity with a 
parabolic temperature profile on the free liquid surface 
with^M* = 350 and v c =» 0. For silicon with (T 0 -T m ) * 
0.05 C and a = 0.5 cm. 
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Figure 4. Computed dimensionless streamlines ij> for surface tension 
driven flow in a floating zone at zero gravity with a 
parabolic temperature profile on the free liquid surface 
with M' » 3500 and v c = 0. Silicon with (T 0 - T ra ) - 0.5 C 
and a * 0.5 cm. 
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Figure 5. Computed dimensionless streamlines for surface tension 
driven flow in a floating zone at zero gravity with a 
parabolic temperature profile on the free liquid surfac 
with M' * 7000 and v c ■* 0. Silicon with (T 0 -T m ) “1.0 
and a ■ 0.5 cm. 
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D. Fluid Flow Coupled with Heat Transfer 

(Electron Beam Heating) 

Although a parabolic temperature profile at the free liquid 
surface may be a reasonable approximation for radiation heating, only 
a complete heat transfer analysis coupled with flow calculations would 
provide us with exact temperature profiles. Electron beam heating is 
commonly used in floating zone melting and was taken as our heading 
mode. Rather than specify the power input, the circumferential temper- 
ature T q at the center of the zone was specified for convenience in 
analysis. The procedure used to solve the coupled heat-momentum 
transport problem was as follows: 

1) The temperature field for conduction was computed by 
neglecting convective heat transfer. 

2) The surface temperature profile from 1) was utilized to 
calculate the first approximate solution of fluid flow. 

3) The temperature field was recalculated using the flow 
fields from 2) . 

4) Steps 2) and 3) were repeated until the temperature and 
the fluid flow fields no longer changed appreciably. 

The resulting streamlines with M' = 35 are drawn in Figure 6 
in which we took values for silicon with the surroundings at the 
melting point, i.e., the heat shielding about the zone is extremely 
effective. While this does not correspond exactly with reality, it 
does show the correct features. Comparing with Figure 2, we see that 
the vortex centers are shifted nearer to the heat source from the 
liquid/solid interfaces, 'ihis is because the steepest temperature 
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gradient is at the center of the zone. Further increasing M 1 up to 

350. by increasing (T -T ) with other conditions kept constant, 

o m 

resulted in the streamlines and isotherms shown in Figure 7. The 
maximum velocity for silicon was 0.55 cm/sec. Comparing the tempera- 
ture fields with surface tension driven flow for M* = 350 with those 
for pure conduction , there was no significant change except for a 
slight one near the center of the zone, as shown in Figure 8. This 
indicates that conduction is much greater them convective heat trans- 
fer as would be expected for the small Prandtl number for silicon. 
(Prandtl number is a measure of the degree to which convection in- 
fluences heat transfer.) The vorticity fields for M' = 350 as in 
Figure 7 are shown in Figure 9. The maximum vorticity and its loca- 
tion for various conditions are also summarized in Table X. 

E. Influence of Prandtl Number 

The Prandtl number was increased from 0.023 to 0.3 and 2.0 by 

increasing the specific heat, and keeping the viscosity and the thermal 

conductivity of the melt constant. As the Prandtl number increases, 

the temperature gradients along the free liquid surface (for fixed 

T -T ) increase near the heat source and also near the liquid/solid 
o m 

interfaces, but significantly decrease inbetween, as illustrated in 
Figure 10. Since the convective heat transfer becomes more significant 
as a result of increasing the Prandtl number, the isotherms are in- 
creasingly distorted from those of pure conduction. This, in turn, 
causes changes in the flow field. For example, the center of the vortex 







M' 


0 


Figure 


M’ * 350 

(CONDUCTION) (CONVECTION) 


. Comparison of isotherms 0 » (T-t )/(T -T ) without 

m o m 

convection and with convection for M' ■ 350 (see 
Figure 7 for the flow field) . 
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DIMENSIONLESS TEMPERATURE 
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<*-v/ <vv 


. Axial dimensionless temperature profiles at the free 
liquid surface; (I) as solved from the electron beam 
beam heat transfer analysis coupled with the convection 
field for M’ = 350; (II) assumed parabolic. 


Figure 10 
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cell near the liquid/solid interface shifts closer to the interface 
as the Prandtl number increases, as shown in Figures 11 and 12. 

F. Influence of Zone Travel on Hydrodynamics 

The effect of zone motion on the flow field was found to be 
negligible in all of our calculated streamlines for zone travel rates 
of up to 5 cm/hour in silicon. However, the effect becomes significant 
when the zone travel rate becomes comparable with the velocity of 
surf ace-driven flow. 

The flow field shown in Figure 6 for M' = 35 (our weakest flow) 
is taken as a model case in order to show the influence of zone travel 
on hydrodynamics in the floating zone. The maximum velocity of the 
melt for M' = 35 is 0.07 cm/sec (or 250 cm/hour), and the average 
velocity is about 70 cm/hour. In Figures 13 through 16, the stream- 
lines in floating zone melting of silicon at zero g are shown for 

freezing rates, 10, 30, 50, and 70 cm/hour with M* = 35 (T -T = 

o m 

0.005°C for our silicon example) . The lower vortex cell floats away 
from the bottom interface and its size is reduced as the zone travel 
rate increases. Thus, the flow patterns with zone motion are quite 
different from those with no motion when the zone travel rate is com- 
parable with the surface driven flow. 




Figure 11. Computed isotherm 6 and streamlines iJj for surface 

tension driven flow in a floating zone at zero gravity 
from the computed temperature fields for electron beam 
heating, with the same conditions as in Figure 7, 
except Pr = 0.3 
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Figure 12. Computed isotherms 0 and streamlines i p for surface tension 
driven flow in a floating zone at zero gravity from the 
temperature fields computed for electron beam heating with 
the same conditions as in Figure 7 , except Fr > 2. 
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Figure 13. Computed streamlines \p for surface tension driven flow 
in a floating zone at a zone travel rate of 10 esn/hour 
under the conditions of Figure 6. 
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Figure 14. Computed streamlines i|> for surface tension driven flow 
in a floating zone at a zone travel rate of 30 cm/hour 
under the conditions of Figure 6. 
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G. Impurity Segregation 

Steady state impurity concentration fields in the melt were 

calculated by introducing the convection fields for M' = 350 with 

zone motion. The values used in the present computations ares 

3 3 

Sc » 5.0, p = 2.28 g/cm , p. * ** 2.54 g/cm , k = 0.1, a - 0.5 cm, and 
c f o 

-4 2 

D = 7.6 x 10 cm /sec, which are typical for silicon with aluminum 

or indium impurity (2) . The dimensionless impurity concentration 

fields, C = m/m^, at steady state are shown in Figure 17 for a 

freezing rate v of 1 cm/hour (V = 0.0366) and in Figure 18 for 
c c 

v of 5 cm/hour (V = 0.183) . The concentration fields follow some- 
c c 

what the patterns of the convective flow fields, which implies that 
the convective mass transfer is large compared with the molecular dif- 
fusion. This is, of course, to be expected from the relatively high 
Schmidt number for silicon. The dimensionless impurity concentration 
profiles in the axial direction at R = 0, 0.5 and 1.0 are plo ted in 
Figure 19 for freezing rates of 1 cm/hour and 5 cm/hour. The magnitude 
of the mass transfer resistance is indicated by the rate of change of 
concentration with distance. Note that it is higher at the center of 
the zone because of the poor mixing between the top and bottom cells. 
This is to be contrasted with boundary-layer flow in which the mass 
trahsfer resistance is concentrated near solid/liquid interfaces. 


* 

When zone melting is initiated the impurity concentration in the melt 
increases (fcr k Q < 1) until the same average concentration is freez- 
ing out of the melt as is being fed into the melt at the other 
interface. It takes several zone lengths of zone movement with 
convective: mixing in the zone until this steady state is reached. 


i 


i 




Figure 17. Steady state dimensionless impur 
tion fields, C = m/m t , in a sili 
a zone travel rate of v c = 1 cm/ 
0.5 cm, M* - 350, Pr f 4 0.0;£3, Sc 
0.1, and D = 7.6 x 10 cm /sec. 
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The impurity concentration in the growing crystal was calcula- 
ted by multiplying the concentration in the melt at the freezing 
interface by the interfacial distribution coefficient, k Q . The result 
is shown in Figure 20. For k Q «* 0.1, the convection in the zone with 
N* « 350 leads to a 10% radial variation in the impurity concentration 
in the resultant crystal at a freezing rate of 1 cm/hour and a 40% 
radial variation at a freezing rate of 5 cm/hour. Note that 5 cm/hour 
is a typical crystal growth rate for silicon. 


H. Constitutional Supercooling 

The criteria for constitutional supercooling indicates that it 
is favored by a decreased temperature gradient and by an increased 
interfacial concentration gradient. The Ilusselt number Nu corresponds 
to a dimensionless temperature gradient and the Sherwood number Sh 
corresponds to a dimensionless concentration gradient, defined as 
follows: 


Nu x 



(1) 


Sh 


3C 

3Z 


+ £./a 


( 2 ) 


From Figures 7, 17 and 18, Nusselt numbers and Sherwood numbers were 
calculated and are shown in Figure 21. Constitutional supercooling 
at the freezing interface is more likely at the center of the crystal 
than at the periphery since the Sherwood number is higher at the 



DIMENSIONLESS IMPURITY CONCENTRATION IN CRYSTAL 



DIMENSIONLESS RADIAL POSITION , R 

Figure 20. Radial impurity concentration profiles in the crystal 
at steady freezing rates of 0, 1, and 5 cm/hour cal- 
culated from the interfacial concentrations in Figures 
17 and 18 for an inter facial distribution coefficient 

of k -0.1 
o 




NUSSELT NUMBER, Nu SHERWOOD NUMBER, Sh 
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FREEZING MELTING 

INTERFACE INTERFACE 

Figure 21. Radial variation of dimensionless temperature gradient 
(Nusselt number) and impurity concentration gradient 
(Sherwood number) at the melt/solid interfaces calcula- 
ted from Figures 7, 17 and 18. 
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center while the Nusselt number is lower at the center. This tendency 
becomes more pronounced as the freezing rate increases since the radial 
variation of interfacial impurity concentration gradient as well as 
absolute interfacial concentration gradient increases, as shown in 
Figure 21. We have found previously that the preferred location for 
interfacial breakdown is strongly influenced by convection (3,4). 

I. Influence of Gravity on Flow 

In a vertical silicon melt with M' = 350 (T -T = 0.05°C for 

o m 

2 

our silicon example) , we found that at 1 g & (980 cm / sec ) the flow and 
temperature fields in the floating zone do not change appreciably 
from those at zero g. When the acceleration is increased to 10 g^, 
the lower vortex contracts while the upper one is expanded, par- 
ticularly at the center of the zone where the flow field is relatively 
weak. This is shown in Figure 22 and may be compared with Figure 7 
at zero g. This means that surface driven flow predominates even on 
earth for many materials for radiant and electron beam heating. It may 
even be important when induction heating is employed and may account 
for some of the compositional inhomogeneities observed. 

J. Influence of Thermal Environment on Isotherms 

The previous results on streamlines were based on low heat flux 
conditions accomplished by increasing the temperature of the surround- 
ings sufficiently or by effective heat shielding. In addition the 
solid phase was not included in the heat transfer analysis. In order 
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to see how much difference this makes, pure conduction analyses 

were made. A low ambient temperature (i.e., T = 23°C) was used. 

c 

Three different thermal boundary conditions at the end of the ingot 
were also investigated. Three possible conditions were considered; 
that is, radiation, insulation and rapid heat transfer from the end 
to the surroundings. For electron-beam heating, the heat flux varied 
substantially depending on the thermal boundary conditions. The 
required power supply computed from the temperature gradients is 236, 
228 and 492 watts for the case of radiation, insulation and infinite 
heat transfer coefficient at the top surface, respectively. This 
result is found to be close to Donald's (5) experimental result (198 
watts) for electron-beam floating zone melting of 1 cm diameter 
silicon. In spite of the changes in heat flux, the isotherms in the 
zone vary little when the length of the solid ingots is the same as 
that of the zone, as shown in Figures 23 through 25. This was also 
true with a reduced length of solid, as shown in Figure 26. More 
significantly, isotherms in the zone with high heat fluxes are not very 
different from those with very low heat fluxes. 

K. Influence of Materials Parameters 

Inspection of the governing differential equations in Appendix 

A shows that the parameter which determines the vigor of the surface 

2 

driven convection is M' = p0(T o -T m )a/y . The temperature difference 
< V T m ) along the melt surface increases as the melting point increases 
and as the thermal conductivity decreases, as well as with decreasing 



ited 

ente 

1410 ' 

ition 






INSULATION 



Figure 24. Computed isotherms, 0 «= (T-T m )/(T m -T c ) , in the ingot above 
the zone center with the fame conditions as in Figure 23, 
except an insulation boundary condition at the end of the 
solid. T 0 -T m = 190°C. 
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M' = 0 M' = 350 

(CONDUCTION) (CONVECTION) 

Figure 26. Computed isotherms, 0 = (T-T m )/ (T m -T c ) , in the solid and 
the liquid above the zone center with the same conditions 
as in Figure 23. The left half of the figure is for pure 
conduction, and the right half is with the convection 
shown in Figure 3. 





heat shielding. Unfortunately we are unable at this time to express 

this quantitatively, although our conputer programs are capable of 

estimating any particular case. If we remove (T -T )a from M' we 

o m 

2 

are left with the material parameters po/U . Values for selected 
materials are shown in Table II, from which it is seen that our test 
material, silicon, is one of the materials more sensitive to tempera- 
ture variations. The parameter which expresses the sensitivity of 
segregation to flow conditions is the Schmidt number Sc. Note that 
Sc is lowest for silicon among all of the materials in Table II. 

While general conclusions are difficult, it seems clear that surface 
tension driven flow will manifest itself in nearly all materials with 
a free melt surface subjected to a temperature gradient. 



TABLE II. SOME PROPERTIES OF LIQUIDS NEAR MELTING POINT 



700 0.06 3.05 2.72 0.024 0.30 2.2 0.34 40 865 
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III. TRANSIENT FREE CONVECTION AT LOW G 
ABOUT A SUSPENDED CRYSTAL 

(Chong K. Kim) 

Definition of Symbols 

a Radius of sphere, cm 

d Diameter of sphere (= 2a) , cm 

D/Dt Substantial derivative. Rate of change as one 

moves along with fluid 

3 3 3 3 

= ttt + u + u +u t— 

3t x 3x y 3y z 3z 

in cartesian coordinates 

2 

g Acceleration, cm/sec 

2 

g Earth's gravitional acceleration, 980 cm/sec 

e 

Gr Grashof number 

d^ggAW buoyancy force 
^2 ~ inertial force 

K Mass transfer coefficient, cm/sec 

2 

p Pressure of the fluid, dynes/cm 

2 

p' Pressure due to fluid motion, dynes/cm 

2 2 

P Dimensionless pressure, p'a /pv 

r Radial distance from the center of the sphere, cm 

R Dimensionless radial distance, r/a 

Ra Rayleigh number, GrSc 

Sc Schmidt number 

V ■> momentum diffusivity 
D molecular diffusivity 

Sh Sherwood number, dimensionless composition gradient, 

Kd/D ■ |3$/3 r| 

t Time , seconds 
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u Radial velocity; cm/sec 

U Dimensionless radial velocity; ua/v 

v Velocity in angular direction, cm/sec 

V Dimensionless radial velocity, va/v 

x,y,z Coordinates in rectangular system 

w Mass fraction of solute 

a Concentration densification coefficient 

6 Angle from negative of accelerational vector 

H Viscosity of fluid, g/cm sec 

2 

V Kinematic viscosity, y/p, can /-ec 

3 

p Density of the fluid, g/cm 

2 

T Dimensionless time, tv/ a 

<J> Azimuthal angle 

$ Dimensionless concentration, W-W -W^ 

ip Dimensionless stream function (Eq. B-8) 

0 ) Dimensionless vorticity (Eqs. B-8 and B-ll) 

Subscripts 

o Surface condition 

» Initial condition 

e Earth condition 

6 Angular 

i,j,m Arbitrary integer 

r Radial 
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T Constant temperature condition 

$ Azimuthal angular 

A. Introduction 

The object of this project was to elucidate the nature of the 
convection about a suspended crystal growing from a solution in a 
spacecraft environment. The small fluctuating accelerations will 
give rise to transient free convection currents. There have been 
very few studies of steady state free convection at very low Grashof 
numbers, and none have been done on transients. (Grashof number is 
roughlj equivalent to the ratio of buoyancy to viscous forces.) 

B. Model 

For simplicity, our attention was devoted to mass transfer about 
a sphere suspended in a solution. The following assumptions were 
employed : 

a) Two-dimensional flow, i.e., no azimuthal flow. 

b) Constant properties (p, D, p) 

-3 -6 

c) Low gravity (g =10 g g ~ 10 g^) 

d) Influence of rate of change of crystal size was neglected. 

The partial differential equations and boundary conditions for 
flow and for mass transfer are shown in Appendix B. Finite difference 
methods were used on the computer. Although the results are expressed 
in dimensionless form, it is helpful to think in terms of a concrete 

example. He considered a 1 cm diameter KC1 sphere growing from water 

o =4 

at 25 C with a supersaturation of 0.1 mol e/l at g ■ 1.25 x 10 g = 

© 

I 

i 

* 

f 


“ f 
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2 

0.122 cm/sec . The pertinent parameters are d * 1 cm, W W o -W w 

- -0.0075, p - 1.18 g/cm 3 , y s 0.01 g/cm sec, a = (Sp/SwJy/o « 0.616 

and D - 1.8 x 10 ^ cnw'sec. These yield Sc = v/D = 465, Gr ** 

3 9 

|d gotAW/v | = 7.9, and Ra « 3670. The boundary-layer approximation 

4 

is expected to be valid only for Ra > 10 £6). In other words, the 
perturbed region in the solution is expected to be on the order of the 
same size as the growing crystal. 

C. Convection-Free Mass Transfer (Diffusion) 

Before studying mass transfer with convection, it is interest- 
ing to first look at the diffusion problem. We consider a sphere in 
an infinite body of solution at zero g. The pertinent differential 
equation is: 

iiil (3) 

3t V/D I 9r 2 R 3 rI 1 ' 

with boundary conditions; 

at T ■ 0 , <(> ■ 0 

at R * 1 , $ ° 1 

at R * 00 , <p « 0 

The solution to this is: 

. 1 R-l 1 , R-l 

<t » — erfc ■ r erfc — — (4) 

R 

This is plotted in Figure 27. The Sherwood number is 
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Figure 27. Concentration profiles for diffusion without convection 
from sphere of radius a into infinite fluid. Dimension- 
less concentration is 1 - <j> = (V. 1 )/(W -W ). Plotted 
_ _ o » o 

from Eq. (4 
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-2Ka 

D 



( 5 ) 


which is plotted in Figure 28. From these figures it is seen that 

2 

to reach steady state requires Dt/a ~ 5, or for our example 

parameters, t ~ 69,600 sec = 19.3 hour. This is to be compared with 

Siegel's (7) predictions based on transient boundary-layer equations 

for Sc » 1. He estimated that the time to reach steady state with 

1/2 

convection on a vertical flat plate is 5.24 (0.952 + Sc) 

|gaAwJ ' d ' = 4760 sec = 1.32 hours for our example. The time 

over which the above solution for pure diffusion should be valid was 
estimated to be 1.80 (1,5 + Sc) (gotAW) d^^ = 1640 sec = 17 min. 


D. Convective Mass Transfer 

In order to reduce the computation time and to approximate a 
realistic geometry, a spherical boundary at ten times the sphere radius 
was chosen. The fluid velocity and the concentration gradient at the 
boundary were set equal to zero. Unfortunately computations were not 
carried on sufficiently long for the convection to appreciably 
influence the mass transfer. As shown in Figures 29 through 34, the 
calculations did show development of a circular roll cell, similar in 
appearance to a smoke ring. This agrees with the steady state results 
calculated by Suriano, Yang and Donlon (8) for a heated vertical plate 
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Figure 30. Constant concentration curves at T = 1.24. Values 
on curves cure -log <p. 
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RADIAL POSITION, R= r / a 


Figure 32. Vorticity as a function of radial position at 54° 
and 90 from vertical for T = 1.24. 




TANGENTIAL VELOCITY, IO*V 
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Figure 33. Tangential velocity at 0 = 90° as a function of 

radial position for three times. For 1 cm diameter 
sphere and a dilute aqueous solution a value of 
V = 0.01 corresponds to about 0.002 mm/second = 

7.2 mm/hour. 
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Figure 34. Radial velocity at 0 = 90° at three times. 
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with a Prandtl number of 10 and Grashof numbers of 0.5 and 1.0, and 
by Hassain and Gebhart (9) for a heated sphere with a Prandtl number 
near 1 and Grashof numbers of up to 0.7. 


The Prandtl number plays the role in heat transfer that the Schmidt 
number plays in mass transfer. 
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IV. EXPERIMENTAL SIMULATION OF CRYSTALLIZATION IN A 
LOW G SPACECRAFT ENVIRONMENT 

(P. J. Shlichta and P. S. Chen) 

Definition of Symbols 

3 

C c Solute molar concentration in crystal, moles/ cm 
Cp Heat capacity of solution, erg/g°K 

2 

D Solute diffusion coefficient in solution, cm /sec 

2 

g Acceleration, cm/sec 

2 

g g Earth's gravitational acceleration, 980 cm/sec 

Gr^ Grashof number for temperature differences, 

L 3 Op/3T) (T^Tjg/pv 2 

Gr Grashof number for composition differences, 

m 3 2 

L (3p/3W) (W o -WJg/pv 

h Heat transfer coefficient between crystal and bulk 
solution, erg/cm sac K 

Heat of crystallization, erg/g 

k Thermal conductivity of solution, erg/cm sec°K 

K Mass transfer coefficient, cm/sec 

L Width or diameter of crystal, cm 

Molecular weight of solute, g/mole 

2 

(°A^o Mass flux of solute into crystal, g/cm sec 

2 

(n ) Mass flux of solvent into crystal, g/cm sec 

Nu Nusselt number, hL/k. Also, dimensionless 

temperature gradient at crystal surface. 

Pr Prandtl number, vpcyk 

Jr 

s Parameter dependent on (Sc/Pr) , from Ref. (11) . 

Sc Schmidt number, v/D. 
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Sh Sherwood number, KL/D. Also, dimensionless 

concentration gradient at crystal surface. 

o 

T Temperature, K. 

V c Linear crystal growth rate, cm/sec. 

W q Mass fraction solute at surface of crystal. 

Mass fraction solute initially in bulk solution. 

a (9p/3w)/p 

3 (3p/8T)/p 

2 

v Kinematic viscosity, cm sec. 

p Density of solution, g/cm 3 

p Density of crystal, g/cm 3 

c 

Ap Change in density due to composition change, 

W g/cm 3 

Ap T Change in density due to temperature change, 
g/cm^ 

A. Purpose and Scope 

The objectives of this portion of the program are a) verifica- 
tion of the theoretical predictions of Section III by observation of 
the growth (or dissolution) rate and convection around suspended 
crystals or polycrystalline spheres at low Grashof numbers, and b) 
observation of crystal growth phenomena at low Grashof numbers for the 
purpose of predicting the feasibility, advantages, and problems of 
growing crystals in a spacecraft environment, with particular emphasis 
on the effect of convection and/or accelerational transients on the 
formation of crystal defects. 
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These goals are to be accomplished by the study of so-called 
low-gravity-simulation (LGS) systems — solute-solvent pairs in which 
the heat of crystallization and the change of density with concentra- 
tion are both small. We had previously demonstrated that for a given 
crystal size and growth rate, LGS systems crystallize at 1 g much as 
normal solution systems would crystallize in a spacecraft environment 
(10) . These systems are being studied by the use of schlieren and 
other optical techniques to observe the steady-state and transient 
convection patterns around growing or dissolving crystals and the 
change of size (or weight) of these crystals as a function of time. 

Some attention may also be given to the correlation between the con- 
vective history and the location of defects in the grown crystal. 

The objectives may be restated as the following specific tasks. 

1) Search for, verification, and characterization of LGS 
systems. 

2) Development and refinement of optical apparatus and 
experimental techniques. 

3) Observation of the dissolution and growth of polycrystal- 
line spheres, immersed in a uniform static solution, by measurement of 
size and/or weight change versus time, time lag for onset of con- 
vection, and qualitative steady-state convection pattern. 

4) Extension of 3) to include mapping of flow velocity and 
concentration gradient patterns as a function of time. 

5) Repetition of 3) and 4) using single-crystal spheres, with 
special attention to the flow perturbations arising as the growing 
crystal changes from a spherical to a faceted shaped 
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6) Experiments similar to 3) and 4) in which the suspended 
sphere is replaced by a large flat plate. 

7) Qualitative observations of the flow pattern around a 
suspended growing faceted crystal with special attention to the pos- 
sible development of spontaneous transients and their correlation with 
the development of crystal defects. 

8) Extension of some or all of the above to include Grashof- 
number transients, induced by accelerations! o* - forced-convection 
pulses. 

9) Low Orashof-number crystal growth under steady-state forced 
convection. 

It is expected that during the present contract activities 
will be confined to tasks 1) , 2) , 3) , 4) and possibly 71 . 

B, Low-Gravi y-Simulation Systems 

The principal effect of gravity upon crystallization is the 
development of free convection currents around the growing crystal. 

The process of crystallization can change the density of the sur- 
rounding fluid in two ways: 

a) by removing solute from the fluid, thereby chang ii g the 
latter' ••• composition, and 

b) by absorbing or emitting heat into the flm.id, thereby 
changing its temperature. 

For either or both of these reasons, the density of the surrounding 
fluid is, in almost all cases, reduced and a rising convecv.cn cur- 
rent is established around the growing crystal. 
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The most commonly used index of the degree of free convection 
in a system is the Grashof number, which is the ratio between the 
buoyant and viscous forces in a system. If we consider only com- 
position difference-, the Grashof number can be defined ass 



( 6 ) 


where 


L = the "characteristic length" of the system, here 
assumed to be the diameter of the sphere or 
crystal. 

3p, - the change of solution density with concentration 

(solute mass fraction) in the vicinity of 
saturation. 

W = the solution concentration (solute mass fraction) 
at the surface of the crystal, assumed to be 
equal to the saturation concentration. 

W„ = the bulk solution concentration, 
g - gravitational (or centrifugal) acceleration 
p = the density of the solution 

V = the kinematic viscosity (viscosity/density) of 
the solution. 



It is obvious that, all quantities except g being held constant, the 
transfer of a crystal 1 izing system from the earth's surface to a 
spacecraft environment has the effect of reducing the system's Grashof 
number by several orders of magnitude. It is also evident that 
alternatively one might reduce the Grashof number — and therefore, 
hopefully, simulate spacecraft-environment crystallization — by 


V 
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a) reducing the size of the system (L) , 

b) reducing the supersaturation (W -W ) and, therefore, the 

oo o 

growth rate, 

c) increasing the viscosity (V) , or 

d) decreasing 3p/9W. 

From a purely fluid-mechanics viewpoint, all of these changes would 
be equivalent. Unfortunately, the process of crystallization involves 
several size-oependent parameters, such as growth step height and the 
ratio of surface energy to bulk energy, so that microscopic crystal- 
lization at 1 g g is probably not a valid simulation of macroscopic 
low-gravity crystallization. Crystallization also involves supersatura- 
tion-dependent phenomena, such as growth step nucleation and surface 
diffusion, so that the extremely slow growth of crystals at very low 
supersaturations is a questionable simulation of faster growth at 
lew-gravity. Finally, order -of-magnitude increases in viscosity cannot 
be effected without substantial amounts of macromolecular additives 
which, as the sad history of attempts to grow crystals from gels 
attests, create problems of their own. Therefore, it would appear that 
most reliable simulation of low-gravity crystallization would be 
attained by using solution systems in which the change of density with 
solute concentration approaches zero. 

We have thus far neglected temperature differences, for which 
the Grashof number is defined as 



( 7 ) 
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in which 9p/2T is the change of solution density with temperature and 
T and T are the temperatures at the growth interface and in the 

O 00 

bulk solution. We can estimate the relative driving forces of heat- 

and mass-transfer by examining very slow transport from a planar 

surface. At steady state the rate of latent heat liberation, 

V c P c AH c , must be equal to the heat flux back into the solution, 

h(T -T ) , where V is the linear crystal growth rate, p is the 
o 00 c c 

density of the crystal, and h = kNu/L is the heat transfer coefficient. 
The growth rate in turn is given by 


W -W 


VC = K — 
c c 


M A l-W 0 [(n A+ n B )/n A ] 


( 8 ) 


where C c is the molar solute concentration in the crystal, K = DSh/L 
is the mass transfer coefficient without interfacial flow, p^ is the 
density of the fluid, M A is the molecular weight of the solute, and 
(n A > o and (rig^ are the mass fluxes of solute and solvent into the 
crystal from the solution. For a solvent-free crystal, = P C / M A 
and (n ) =0. Inserting Eq. (8) into the steady state heat balance 

we obtain for a solvent-free crystal 


h (T -TJ = V p AH = Kp AH (W -W )/(l-W ) (9) 

o 00 cc c f c 00 o o 


The ratio of the density change due to the compositional change to 
that due to the temperature change is therefore, 
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Ap„ a(W -W ) Gr h(l-W ) 

H W o °° m a o 

Ap T = B(T o -TJ " Gr h = 6 kp f AH c 
From Ref. (11) we note that 

Sh = Nu/s s Nu/0.737 (Pr/Sc) 1/3 , 
and this with the definition of Sh and Nu gives 


h k Sc ln /Pr\ 

— = — ~ Pl C — (0.737) ( — ) 
K D K f p Pr \Sc/ 


= 0.737 p £ C (Sc/Pr) 
f P 


2/3 


Substituting this into Eq. (10) we obtain 


( 10 ) 


(ID 


( 12 ) 


Gr C /c . \?/3 

gT- 0 - 737 «-V < l3 > 

h c ' 

Taking c p = 1 cal/g°C and (Sc/Pr) = 140, we obtain for (CH 3 ) 4 NC1/H 2 0 
(a typical LGS system with a AH of 2.4 cal/g) Gr m /Gr h = 270, showing 
that the density change due to thermal effects is negligible. How- 
ever, the ratio of the driving forces due to compositional and thermal 
effects is, from Eq. (13) and Ref. (11), 

SGr /„ \l/3 C 

<0.544, jg-U-W 0 ) <14, 

fl ' c 


a 


38 for (CH 3 ) 4 NC1/H 2 0 



1 


72 


Therefore, although the Grashof number for the temperature difference 
is estimated to be much less than the Grashof number for the compo- 
sition difference, it might have some adverse effect on the simulation 
of extremely low-gravity crystallizations. It is therefore desirable 
to choose LGS systems for which AH c is as small as possible. Also it 
is obvious that the simulation of growth from melt at low-gravity is 
unfeasible, since the Ah of melting is always considerable. 

Preliminary experiments, reported elsewhere (10) , disclosed 
that the tetramethylammonium chloride/water system, because of its 
extremely low 3p/3w, low AH, and high viscosity, appeared to be a 
valid LGS system. Unfortunately, its extreme hygroscopicity and 
tendency toward dendritic growth made further research with this system 
unattractive. 

Accordingly, at the beginning of the present contract, a 
literature search was initiated to find other LGS systems, using the 
following criteria: 

1. The change of solution density with increasing solute con- 
centration, in the vicinity of saturation, must be as close as possible 
to zero. This is equivalent to requiring that the density of the 
solvent equal the density of the dissolved solute and is approximately 
equivalent to requiring that the solvent and solid solute have the 
same density. Accordingly, the latter requirement was used in search- 
ing for candidate solute-solvent pairs. 

2. The system must have a low heat of crystallization. It is 
still not certain how strict this requirement must be. 
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3. The system must yield good equidimensional single crystals. 
(Our early fears — that all LGS systems would tend toward dendritic 
crystallization — were rapidly dispelled.) 

4. The solute must melt well above room temperature and have 
a low vapor pressure. 

5. The solvent must boil far above room temperature and have 
the lowest possible vapor pressure, so as to minimize the likelihood 
of evaporative nucleation at the surface of the solution. 

6. Both solvent and solute must be stable and substantia'.ly 
nontoxic with regard to prolonged exposure to vapor and occasional 
skin contact. 

7» The room-temperature solubility should be greater than 
about 0.05 weight fraction, to permit reasonable growth rates, and 
lass than about 0.20 to facilitate comparison with theory. Also, 
solubility should increase with temperature , so as to facilitate 
the preparation of supersaturated solutions. 

8. Ideally, the solvent should be a variable mixture of two 
chemically similar compounds of slightly different density, such as 
H^O-D^O. This permits the variation of 3p/3w, so that the Grashof 
number of the system can be "tuned" to zero or very low values. 

9. The system should possess one or more of the following 
characteristics to facilitate measurement of the convection and/or 
solute concentration around the growing crystal: 

a) A strong variation of index of refraction with 

concentration so as to permit observation by schlieren 


methods 



c 


74 


b) A "window" in the visible region wherein the 
solvent is transparent and the solute moderately absorbent 
(a - 0.1 cm so as to permit absorption spectrophotometry. 

c) Electrical conductivity and electrolyzability, 
so that convection can be visualized by the generation of 
minute hydrogen bubbles (12). 

10. Both solvent and solute should be readily available and 
reasonably inexpensive. 

A previous search — which eventually uncovered the tetra- 
methylimmonium chloride/water system — was believed to have exhausted 
the like3.ihood of finding a water-soluble inorganic salt with a 
density near 1.0. Therefore, the present search was confined to 
organic solutes and solvents. This was also advantageous since the 
heat of crystallization is usually small in organic systems. There is, 
unfortunately, no published tabulation of organic compounds according 
to density. Therefore, our search commenced with the compilation of 
such a tabulation from standard organic compendia (13,14) . This 
compilation (Table III) is by no means complete or systematic, since 
on the one hand we first considered standard solvents and solutes 
which were known to crystallize well, and on the other hand we rejected 
any compound which failed to satisfy criteria 3) , 4) or 5) by having 
an unsuitable melting point, boiling point, or crystal habit. (Many 
solutes, for example, were rejected because they were known to crystal- 
lize as needles.) This phase of the search is by no means complete 
and it will eventually be resumed. 



TABLE III. SUITABLE SOLVENTS AND SOLUTES TABULATED ACCORDING TO DENSITY 
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(continued) 



TABLE III. (Continued) 


i 



1-bromobutaxnl 



TABLE III. (Continued) 
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Table III was then Inspected for possible solute-solvent pairs 
which were tentatively accepted or rejected on the basis of toxicity, 
availability and solubility. When solubility data was unavailable, 
it was estimated on the basis of chemical similarity or analogy to 
known solubilities. Table IV lists some representative rejected 
systems together with the reasons for their rejection. Table V lists 
the systems which, after the first round of screening, appeared to be 
worth further study. Three of these candidate systems wer.. then 
checked for solubi]ity and for variation of density with concentration. 

A Westphal balance was used for the latter measurement. Known 

weights of solute and solvents were combined in varying proportions to 
provide a set of solutions of different concentrations. Solubility 
was estimated by adding small known weights of solute until apparent 
saturation was reached. 

Density measurements for the first three candidate sys.Jns are 
shown in Figure 35. It will be noted that, of these, the thymol/ 
ethyl carbonate system has a 3p/3w of almost exactly zero. Moreover, 
the dp/dw of the other two systems is low enough to make them also 
potential candidates for low-Grashof -number investigations. 

The crystallization form of these systems was investigated by 
preparing saturated solutions at elevated temperatures (e.g., 50°C) 
and then letting them cool slowly. The camphor crystallized out from 
ethyl carbonate solution as needles. Both thymol systems produced clear 
euhedral crystals up to 2 mm in diameter. 
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TABLE IV. UNACCEPTABLE LGS SYSTEMS AND REASONS FOR REJECTION 
phenols water - toxicity 

citric acidsCCl^ - toxicity and probable low solubility 
borneol: 3-butanoic acid - stench of solvent 
benzoic ncidsglycerin - poor crystal form and low solubility 
terebic acidsbutyl alcohol - unavailability of solute 
acenaphthalenes 0-xylene - poor crystal form (needles) 
acenaphthalene : ethyl carbonate - poor crystal form (needles) 
acrylic acids aniline acetate - very low M.P. and instability 
camphenesM-xylene - ]ow M.P. of solute 


TABLE V. POSSIt T E ACCEPTABLE LGS SYSTEMS 

thymols ethyl carbonate 

thymol : acetylacetone 

camphors ethyl carbonate 

borneol s anisole 

borneol scarbitol acetate 

B-naphthylamines benzyl acetate 

methyl hydroxylamine swater (despite low M.P.) 


I 
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Figure 35. Density versus concentration of 3 candidate LGS systems 
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Therefore, despite their high solubility, which makes com- 
parison with transport theory more difficult, and despite the rather 
high vapor pressure and potential low-level toxicity of the solvents, 
the thymol/diethyl carbonate and 'zhymol/acetylacetone systems appear 
to be potentially suitable LGS systems. Thus far, growth and dis- 
solution have beei studied only in the thymol/diethyl carbonate system. 
The results of these tests will be discussed below. 

C . Preparation and Characterize cion of Solutions 

In order to calculate the Grashof number for a given crystal- 
lization it is necessary to know the density, viscosity, solubility, 
and change of density with concentration for the system at the 
temperature and concentration of crystallization. When the experiment 
involves a solution prior to or during crystallization, or when the 
possibility of convection due to heat transfer ic considered, it is 
also necessary to know the density and solubility as a function of 
temperature. In addition, it is useful to know the index of refraction 
of the solution as a function of concentration and temperature. 
Refractive index measurements provide a rapid means of determining 
solution concentration and an easy method of determining solubility. 
Accordingly, the density, viscosity, and index of refraction of thymol/ 
diethyl carbonate solutions are belrg determined as a function of con- 
centration and temperature. The index of refraction (Figure 36) is 
determined by means of a Bausch and Lomb Abbe refract ometer. The 
refractive index measurements are good to at least three significant 
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figures. It will be noted that there is an appreciable increase in 
index of refraction with concentration, so that the system is shown 
to comply with criterion (9a) . The viscosity (Figure 37) as 
measured by a U-tube capillary viscosimeter, depends strongly on both 
concentration and temperature. Attempts to measure density by usi..g 
a Westphal-balance sinker on a precision microbalance led to un- 
accountably erratic results and these measurements were repeated using 
a pycnometer (Figure 38) . At present, the solubility of this 
system is being determined as a function of temperature by equi- 
librating solutions with excess solute, determining the resultant 
index of refraction, and then extrapolating the data in Figure 36 to 
determine the saturation concentration. In making these measurements, 
as in the preparation of solutions for the first growth experiments, 
a difficulty was noted which is peculiar to LGS systems and which 
would also be encountered in spacecraft experiments. Normally, 
solutions are saturated by first dissolving as much solute as possible 
at a higher temperature and then cooling the solution to the final 
temperature where, in the presence of a few seed crystals of solute, 
excess solute crystallizes out. In ordinary solutions, the process 
of crystallization causes sufficient convection to mix the solution, 
so that storage at constant temperature for a few days to a week, 
with only occasional stirring or shaking, is sufficient to insure 
equilibration. In I/SS systems, however, free convection is almost 
negligible and equilibration can be achieved only by stirring or 
shaking the solution continuously. This precaution is also necessary 
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Wt Fraction Of Thymol In Ethyl Carbonate 


Figure 37. Viscosity of thymol-diethyl carbonate solutions at 21.0°c 
25.7 C, and 30.2°C. 




when determining the heat of solution, a measurement which is now being 
carried out on the thyiuol/diethyl carbonate system. 

D. Preparation of Crystals and Polycrystalline Spheres 

Thus far, the procurement of single crystals has presented no 
problem. Cubic cleavages of KBr and KC1, cleaved from melt-grown 
crystals, appear to be thoroughly satisfactory. NaClO^ crystals up 
to several centimeters on a side may readily be grown from solution. 
Thymol crystallizes so easily from melt or vapor that commercial 
material is often sold as single crystal lumps. Moreover, clear 
euhedral crystals can be grown from the LGS solution. 

Polycrystalline spheres were produced by pressing powder between 
hemispherical steel dies. The original die design (1/2" hemispheres 
cut into 5/8" diameter rods as shown in Figure 39) resulted in equa- 
torial ridges and stress concentrations which often caused the sphere 
to fracture into halves. The dies were then modified by enlarging the 
hemispherical surface to the same diameter as the rod and by cutting 
off the ends of the rods so that when the die rods barely touch, 
forming a spherical cavity, their ends are flush with the confining 
ring (Figure 40) . 

The spheres are prepared by pouring into the die chamber a 
weighed amount of powdered solute, corresponding exactly to the 
theoretical weight of a finished sphere. The upper die is then in- 
serted and pressure gradually applied until an abrupt increase is 
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noted. (This indicates that the die rods are flush with the con- 
fining ring.) The pressure is maintained for about 5 minutes and 
then removed. The die rods are then pushed through the ring and 
pulled apart and the sphere is removed. It is sometimes advan- 
tageous to reposition the sphere between the dies and repeat the 
pressing operation several times; this produces a more symmetrical 
sphere. 

In this manner, translucent coherent spheres of KBr, KC1, 

NH^Cl, and camphor can be readily produced by room temperature press- 
ing. Thymol is too brittle to press well at room temperature, but 
produces a coherent sphere at 45-50°C. Thymol also tends to adhere 
to the dies, so that uhe sphere has to be removed from between the 
die rods by gently heating them. 

Crystals and spheres are usually mounted by drilling a #80 hole 
through them and then passing a fine wire through the hole. During 
growth and dissolution experiments, the wire is suspended horizontally 
across the center of the cell. Thymol crystals and spheres can also 
be mounted by melting them onto an electrically heated wire. In 
future experiments, the wire will sometimes consist of a chromel- 
alumel composite with a junction at the center or surface of the sphere 
and in some cases a second junction at the side of the cell, so that 
thymol gradients can be measured directly. Thus far, the only solute 
which has presented any difficulties is NaClC> 3 , which is too brittle 
to press or drill at room temperature. 



E. Optical and Photographic Apparatu s 


The preceding theoretical sections suggest that the most im- 
portant observations during growth or dissolution would be of a) the 
average growth rate, as measured by the change of size or weight 
of the crystal as a function of time, b) the qualitative form of the 
convection pattern around the crystal, and c) the convective velocity, 
especially as a function of time. Similarly, the second objective 
of the experimental program could best be met by continuous observa- 
tion of the convection pattern around a growing crystal. These 
requirements suggested time lapse schlieren photography as the tech- 
nique of choice. Therefore most of our activities to data have been 
concerned with the assembly, testing and refinement of a 4-inch concave 
mirror Topler schlieren apparatus with attached time-lapse movie camera 
The Topler schlieren method provides an image in which the 
brightness at each point corresponds to the magnitude and sign of the 
horizontal (or vertical) component of the index -of-refraction gradient. 
The apparatus is shown schematically in Figure 41. Light from a point 
source (A) is collected by a lens (B) and formed into a parallel col- 
limated beam which passes through the sample cell (C) and through a 
second lens (D) which refocuses the undeviated light at point E. Here 
a knife edge cuts off part of the focused beam, thereby uniformly 
dimming the background of the image. Horizontal refraction gradients 
deviate the light to the left or right of the main focal point. 
Accordingly, these parts of the image appear either darker or lighter 
than the background. Thus, rising convection columns of lower density 
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(and hence less refractive) fluid are revealed on the schlieren image 
by an apparent (but deceptive) three-dimensional relief. So it is 
seen that, in addition to providing a semiquantitative measurement of 
dn/dx, the schlieren technique is also a means of flow visualization 
and can be used to measure eddy velocities and similar phenomena. 

However, it does not provide any measurement of flow velocity, although 
it is fairly easy to incorporate a flow marking device into a schlieren 
apparatus. Moreover, a grid of opaque lines or dots may be introduced 
between B and C; the deviations of these grid points on the schlieren 
image may be used to provide an additional measurement of the index- 
of-refraction gradient. Final?/, with the knife edge removed altogether, 
a sharp silhouette image results j this "shadowgraph" mode is very useful 
for measuring changes of size in tha sphere or crystal. 

In practice, in order to minimize chromatic and spherical 
aberation, concave mirrors are often used in place of lenses. A typical 
configuration of £his type is s^own in Figures 42 and 43. Ideally, to 
avoid astigmatism, B and D should be off-axis parabolic mirrors; these 
are unfortunately prohibitively expensive. The experiments described 
below were performed with an inexpensive (~ $100) pair of 4-inch 
diameter 36- inch focal length telescope mirrors whose principal virtue 
was their availability. 

The point light source (A) is a Sylvania concentrated-arc lamp 
(Model C-10) in which a tiny (0.4 mm mean diameter) brilliant (47 
candles/mm) arc is maintained between a pointed rod and a concentric 
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Figure 43. Alternative schlieren system using mirrors 
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ring of sheet zirconium. Although the arc can theory wander around 
the circle between the point and the ring, this has thus far caused 
no apparent difficulty. In all the crystallization experiments per- 
formed thus far, the lamp was used as the light source without any 
collimation or spatial filtering. 

The chief disadvantage of this apparatus, as originally con- 
structed, was its poor image resolution which, in the shadowgraph mode 
(with the knife edge removed) was no better than 0.5 mm; in the 
schlieren mode, it was somewhat worse (Figure 44) . On the assumption 
that this was probably due to stray light from the arc lamp, a collimator 
was constructed whereby the arc image was projected, by means of two 2" 
diameter 8" focal length achromatic lenses, onto a pinhole positioned 
at point A. Using a homemade 0.6 mm diameter pinhole, the resolution 
was greatly improved (Figure 45) without substantial loss of image 
brightness. Although now adequate *or our present needs, the . 
tion can presumably be improved further by using an even smaller pinhole. 
Arrangements have also been made to perform some experiments using a 
focused laser in place of an arc lamp; this will involve the use of a 
4" achromatic lens schlieren apparatus which is now being constructed. 

Another serious problem has been optical distortion from the 
solution cells. The plexiglas cells used in the previous experiments 
were shown to have numerous schlieren striae due to minute index-of- 
refraction variations in the superficially perfect plexiglas sheets 
used for their construction. Moreover, when these cells were filled 
with solution and inserted into the schlieren system, they shifted 
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the focal point, indicating that the cell windows had a lens-like 
distortion. 

Two sources of "schlieren free" plate glass were located and 
all-glass cells were constructed by joining edge-ground glass plates 
with epoxy cement. These cells, when inserted into the schlieren 
system, still show a few isolated schlieren striae (presumably due to 
compos i tiona 1 inhomogeneities) but are adequate for present use. 
Moreover, they cause no shift in the focal point of the system. 
Inquiries of vendors indicate that completely schlieren-free cells 
would cost several hundred dollars apiece, so that no further pro- 
curement is contemplated during the present contract. The epoxy 
cement used (Devcon "Two Ton") , although rated best for adhesion to 
glass , was slowly attacked and destroyed by the thymol/ 
diethyl carbonate solutions. A series of trial experiments, using 
glued microscope slides soaked for days in various solvents, showed 
that Dow Corning Silicone Adhesive is the most adhesive and stable 
cement for our purpose. 

Three modes of photographic data recording have been used 
(Figure 46) s 

a) Ordinary (lens focused) photography outside of or trans- 
verse to the schlieren apparatus (e.g., through the side walls of 
the solution cells) . Various Polaroid or 35 mm cameras are available 
for this purpose. Our Bolex time-lapse movie camera has not so feu: 
been used for this purpose. 
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b) Schlieren and/or shadowgraph images can be taken with 
a lensless camera at any point beyond E. Since the magnification 
is a function of distance (M = (E-F) "/36") , a moveable track has 
been constructed so that either overall views or closeups can be 
photographed. Two camera backs, a 35 mm Minoleta ST-102 (chosen for 
its ability to take precisely registered double exposures) and a 
Bolex 16 mm movie camera with automatic time-lapse, are inter- 
changeable in this position. 

c) Shadowgraphs at 1:1 magnification can be taken by in- 
serting a camera back between C and D. Plans are now underway 

to record the dissolution of polycrystalline spheres using the 35 mm 
camera back in this mode while (between still photographs) the 
movie camera is used in mode b) . 

Other, non-photographic techniques are also available. For 
example , a strain-gauge microbalance is available for continuously 
recording the weight of the growing or dissolving crystal or sphere. 
This apparatus may be used in experiments on normal systems, such as 
KBr/water, but it is of very little value with LGS systems, in which 
the imnersed crystals are practically weightless due to buoyancy. 
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Experimental Results 

Thus far, our principal concern has been with tasks (1) and 

| (2), i.e., searching for and evaluating candidate LGS systems and 

I 

| developing necessary apparatus and techniques. Therefore, growth and 

f -- 

r ' dissolution experiments were performed primarily to test solution 

L 

* 

j systems or apparatus rather than to obtain interpretable data. None- 

theless, preliminary experiments, which are summarized in Table VI, 

• have yielded interesting results. 

D issolution of NH 4 CI Spheres : At the beginning of the pro- 

gram, it was observed that a polycrystalline NH^Cl sphere, immersed 
in pure water, dissolved rapidly while the spherical form changed to 
a vertically elongated "nose cone" shape (Figure 47) . Later photo- 
graphs taken during the dissolution of an NH^Cl sphere in a 90% 
saturated aqueous solution showed that the sphere gradually changed 
to a vertically elongated spheroid with an equatorial ridge (Figure 
48) . (The flat bottom is assumed to be due to the breaking away of 
a cracked portion of the sphere.) This observation is qualitatively 
consistent with theoretical expectations based on the assumption of 
downward flow, since the top and bottom of the sphere would be 
stagnation points while the maxima in both flow velocity and concen- 

I tration gradient would be in the upper quadrant. We plan to repeat 

S these experiments over a wide range of Grashof numbers with simul- 

1 

| taneous photographs of profile and schlieren. 

Dissolution and Growth of KBr ; In experiment 40331A, a solu- 
tion of KBr in water was saturated at 30°C, poured into a heated 
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TABLE VI. (Continued) 






Figure 47. 


Successive stages of dissolution of NH 4 C1 sphere 
in water. (Experiment 40102) 



Figure 48. Successive stages of dissolution of NH.Cl 
sphere in 90% solution. (Experiment 40312) 
From putside in, 0 min, 10 min, 50 min, 240 
min, 270 min. 
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solution cell, and cooled to room temperature (20 C) , thereby pro- 
ducing a supersaturation of about 8%. Rapid skeletal growth was 
observed at the corners of a seed crystal cube. Abundant nuclei 
were formed on the cell walls and on the cleaved crystal surface. 
Well-defined laminar schlieren were observed around and above the 
seed crystal. The rate of upward convection of the plume was 
estimated at 1 to 2 mn^sec. After 2 hours, the bottom of the cell 
was filled with tiny acicular crystals and the original crystal 
surface was overgrown with needles and/or dendrites at its edges and 
corners. In a later experiment (40513A) , a solution was saturated 
at about 4° above room temperature, heated to destroy stray nuclei, 
and allowed to cool slowly. A cleaved KBr seed was inserted about 5° 
above room temperature and observed to dissolve rather rapidly so that 
the crystal became asymmetrically wedge shaped (the top being slightly 
thinner than the bottom, as would ne expected) . The crystal was re- 
moved and reimmersed when the solution was 2° cooler. Schlieren 
observation then showed that there was almost no convection around 
the crystal except for slight slowly moving oblique schlieren, pre- 
sumably due to thermal gradients. As the solution cooled further, 
upward convection rapidly developed which, although apparently uniform, 
became violently unstable for several seconds whenever the system was 
disturbed (e.g., by a knock on the supporting table) or the crystal 
displaced. Another experiment (40515) , in which a single-crystal KBr 
cylinder was immersed in a 98% saturated solution, showed an apparently 
uniform downward flow with plume velocities of about 5 to 10 mm per 
second, with similarly violent instabilities whenever the support table 
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was tapped with a finger (Figure 49) . Slight flow irregularities 
were observed to occur occasionally even under apparently vibration- 
free conditions (Figure 50) . During these experiments some slight 
residual turbulent convection, apparently due to the recirculation of 
the main convection current after it encountered the cell walls, was 
also noted but it was far weaker them the primary current around the 
crystal. 

Thymol/Diethyl Carbonate System : In our first experiment 

(40331B) , a solution of thymol in diethyl carbonate was equilibrated 
for 3 days at 30°C, poured into a preheated glass cell, and allowed to 
cool 5°. Insertion of a 1 cm thymol crystal initiated profuse quasi- 
dendritic (i.e., bladed) growth, both on the seed crystal and through- 
out the solution. After 30 minutes of growth, during which a poly- 
crystalline cluster had formed around the seed crystal, a stationary 
vertical convection column was observed extending upward from the 
cluster. This schlieren pattern did not appear to change during 5 
minutes of continuous observation. 

After a day of growth, during which all spontaneously initiated 
crystals had slowly settled to the bottom of the cell, it was noted 
that the later stages of crystallization (i.e., at lower supersatura- 
tion) had produced transparent euhedral crystals. A dissection of the 
cluster similarly disclosed a core of radial blades and an outer region 
of transparent euhedral crystals. 

A subsequent experiment (40513B) under similar conditions was 
monitored overnight by time-lapse cinematography. The next morning a 





Figure 49. Flow about KBr , perturbed by a tap on the 
container. Frames from movie taken at one 
frame per second, s 1 irting at upper left. 
(Experiment 40515) 








large radiating dendritic cluster was found in the cell. Unfor- 
tunately, most of the movie was badly overexposed, due to a shutter 
failure. 

These experiments convinced us that our solutions had too high 
a supersaturation and that our equilibration procedures wen in- 
adequate. We then realized the need for continuous stirring during 
equilibration and modified our procedures, as was described in Section 
C. Subsequent experiments showed that clear euhedral thymol crystals 
could be grown at rates of several millimeters per day from solutions 
cooled 2° or 3° below their saturation point. 

The convection pattern during slow growth or dissolution is 
quite remarkable and is very different from the theoretical predictions 
in Section III. For example, immersion of a polycrystalline thymol 
sphere into a 98% saturated solution (experiment 40514A) caused a 
downward laminar plume to flow at a rate of about 0.3 mm/sec (Figure 
51) . This convection stream was unaffected by disturbance of the cell 
or support table. Even displacement of the crystal merely caused a 
smooth curvature of the convection stream. What is most remarkable, 
however, is that this convection stream remained coherent and sharply 
defined even after it encountered the cell wall which deflected it 
upward. The ultimate effect, a tangle of twisted schlieren, is 
reminiscent of the evolution of a pyrotechnic "snake" (Figure 52) . 

G. Conclusions and Future Plans 

From the foregoing experiments it appears that the apparatus 
and techniques veloped thus far, although capable of considerable 




Figure 51 


Plume from thymol dissolving in diethyl carbonate 
shortly after initiation of dissolution (experiment 
40514). 5 seconds between frames. 




Figure 52. 


Meandering plume from thymol dissolving as ethyl 
carbonate several hours after initiation of dis- 
solution. 'Experiment 40514.) 5 seconds between 

frames. 
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further improvement, are at present adequate for the observation of 
a great deal of interesting phenomena. Therefore, during the next 
six months we plan to give primary attention to the performance of 
experiments and carry out only those apparatus improvements which 
are readily feasible or particularly needed, such as improved shadow- 
graph resolution and cell-temperature control. 

It also appears that the thymol/diethyl carbonate system has 

been shown to crystallize at low Grashof numbers and that, at least 

at reasonably low growth rate, heat effects are not important. This 

is perhaps sufficiently demonstrated by comparing the dissolution of 

a KBr cylinder in 98% saturated aqueous solution with the dissolution 

of a thymol srhere in 98% saturated diethyl carbonate solution. As 

Table vi shows, the convective velocity varies approximately as the 

1/4 power of the Grashof number. Therefore, it seems safe to assume 

that at constant gravity and acceleration, the thymol/diethyl carbonate 

system is a valid model at 1 g for the crystallization of ordinary 

e 

-4 

solution systems in fields of -10 g g . 

There is a discrepancy between the vortex-like cell and the 
convection pattern predicted in Section III and the experimentally 
observed plume "snake". The resolution of this discrepancy will be 
one of the primary goals during subsequent months. At present, we 
tend to regard the reluctance of the convective striae to mix with 
the bulk solution and the re-emergence of these striae, reflected by the 
cell walls, into the region of crystallization, as significant warnings 
of problems that will be encountered during crystallization in space. 

For this reason, we ' .-lieve that the present program should eventually 


T 



give considerable attention to the theoretical analysis and exper- 
imental simulation of low-gravity crystallization with forced con- 
vection. 

The above experiments also demonstrate that our "low- 
gravity-simulation" systems, which might more properly be called 
"acceleration-insensitive" systems, cannot easily be used for pre- 
dicting the effect of accelerational transients on crystallization 
in a spacecraft environment. However, the convective chaos effected 
by accelerational transients on ordinary systems, such as KBr/water, 
has been shown to be so spectacular that it certainly merits further 
study. 

Finally, successful as the thymol/diethyl carbonate system 
appears to be, it seems necessary to keep searching for new systems — 
not only because of the potential toxicity of the present system, 
but also to demonstrate that the effects observed cure not peculiar to 
one system, but are characteristic of any system crystallizing at a 


low Grashof number. 
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V. FUTURE PLANS 

The low gravity simulation experiments will be continued as 
outlined in Section IV. The effort up until now has concentrated 
primarily on experimental methods and techniques. In the future 
emphasis will shift more to actual growth and dissolution experiments. 
Flow phenomena and its influence on crystal growth will be studied. 

If additional funding is forthcoming, a variety of theoretical 
calculations can be performed. Some possibilities related to surface 
driven flow are: 

1. Extension of previous work on floating zone melting to 
other values of the parameters, i.e., higher Pr and Sc, higher 
surface tension driving forces, different equilibrium distribution 
coefficients, etc. 

2. Transient analysis of floating zone melting of high Prandtl 
number material with radiant heating, for which oscillatory convection 
is expected. 

3. Determination of shape of floating zone due to a surface 
tension gradient in zero g. 

4. Surface driven flow and melt shape in the gradient freeze 
method with a free melt surface, as observed in several Skylab experi- 
ments, but not yet explained or understood. 

At best, only two of these tasks could be expected to be completed in 
one year by an experienced Ph.D. research assistant. 

Possible future calculations for free convection about a sus- 
pended object might investigate: 
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1. The influence of the size of the crystal gorwth container. 

2. Influence of Gr and Sc. 

3. Influence of finite solubility (changing diameter of spbar >. 

4. Other geometries, such as cylinder or cube. 

5. Perturbations of a steady state by g fluctuations. 

A graduate student research assistant could be expected to complete two 
of these in a one-year period. 
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APPENDIX A 

ANALYSIS OF SURFACE DRIVEN FLOW IN 
FLOATING ZONE MELTING 


(Chong E. Chang) 

The definitions of symbols and the assumptions were given in 
Section II. The coordinate system and geometry for floating zone 
melting were shown in Figure 1. Following are the fluid flow, heat 
transfer, and mass transfer equations used in the analyses (A-l, 
A-2, A-3) . 


A. Dimensional Equations 

(1) Continuity Equation 

At constant melt density, conservation of mass requires 

that 

3v v 3v 


+ JL + — £ = o 

3r r 3z 


(A-l) 


where v is the radial velocity and v is the axial velocity at r and z. 
r z 


(2) Momentum Equations 

Force balances yield the Navier-Stokes equations for a 
Newtonian fluid: 

/ 3v r 3v r \ 3P (*\ 1 3v r v r ^rl 

Plr3r zTT 3r M ^2 r 3r 2 3j 2 


3v 


3v 


p « V r 37" + V z 37 


Vv . 3v 3‘v 

- tr + M— 4 + + 9 


3P 

dz 


3r 


3z 


(A-3) 



where the total pressure P may be decomposed into a hydrostatic pox .on 


P Q and a dynamic portion p, representing the change in local pressure 
due to the convection. Thus we define 


where 


and 


P = P o + p 


p Q = P f g(*-z) , 


p[l+$(T-T m )] , 


0 si /V| 

P \ 9T /n 
\ / P i 


(A-4) 

(A-5) 

(A-6) 

(A-7) 


m 


The boundary conditions for the flow field are: 


a) At z = £ 


c) At Z = 0 


V » 0, 
r 


V = -V 

z cl p, 



3v 

_ i 

T z 3z 
3v 



d) At r = 0 v = 3* — = 0 
r or 


(A-8) 


e) At r - a 


f) 


v ■ 0 and 
r 


3v 


2 „ I lx * £ 21 

3V li 3z y 3z 


where p is the density of the liquid, p c is the density of the solid, 
p^ is the density of the liquid at the melting point, y is the 
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viscosity, g is the acceleration field opposite to the z-axis, A is 

one half of the zone length, a is the radius of the zone, T is the 

temperature, is the melting point, 3 is the thermal densification 

coefficient, v is the zone travel rate, y is the surface tension, and 
c 

a = 3y/3T is the temperature dependence of the surface tension. 
Boundary conditions (a) and (b) arise from the zone motion, (c) and 
(d) from symmetry, (e) fron confinement of flow to the zone. Bomdary 
condition (e) results from equating the surface tension gradient force 
to the viscous force at the melt surface. 


(3) Heat Transfer Equations 

A heat balance over a differential element yields for 
constant properties: 


P c. 


/ 3t 3t\ ./s t l 3 t 3 2 t\ 

V 3z Vr 3r / = V r 3r 3r 2 / 


(A-9) 


The boundary conditions for heat transfer are: 


3t 


a) At r = 0 , 


b) At r ■ a , 


. 3t , 4 4 

-k r- = C£ (T -T ) 
3r s c 


c) At 


( r = a ) 

< and > , T 

(z - 0) 


d) At z ■ i , T » T 


m 


e) At z ■ -1 , T ■ T or 

m 


3T 


f) At z ■ 0 , ^ ■ 0 (only for v c ■ 0) 


(A-10) 
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where T Is the temperature of the surroundings, T is the melting point, 
c m 

T is the temperature of line heat source at the periphery of the center 
o 

of the zone, c is the heat capacity of the melt, k is its thermal con- 
P 

ductivity, a is the Stefan-Boltzmann constant, and e is the emissivity 

s 

of the melt. Boundary conditions (a) and (f) arise from symmetry, (b) 
by equating condition to the melt surface to radiation into the sur- 
roundings, and (c) , (d) and (e) by definition. 

When the solid phase was included in the heat transfer analysis, 
the following alternative boundary conditions were used instead of the 
boundary condition 10 d) . Radiation, insulation, and infinitely rapid 
heat transfer to the surroundings at the top of the solid were all con- 
sidered. The boundary conditions for these possibilities are: 


d') At z = b , 



4 4 

ae (T -T ) for radiation 
s c 


3t 

3z 


0 for insulation 


(A-10) 


T ■ T for infinite heat 

c transfer coefficient 


where b is the distance from the center of the zone to the end of the 
solid rod. 


(4) Mass Transfer 

An impurity balance over a differential element yields: 


3m . 3m 
v — + v — — 
z 3z r 3r 



(A-ll) 


Boundary conditions for mass transfer in the mass-centered 
reference system are depicted in Figure A-l, and are as follows (A-4) 
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a) At r ■ 0 , 


3m A 

5r " 0 





(A-12) 


where m is the mass concentration of impurity, D is the binary molecular 
diffusion coefficient, m .is the impurity concentration in the crystal, 

CD 

m ct is the inpurity concentration in the feed ingot (assumed constant) , 
k is the impurity concentration in the melt at the growing inter- 
face, m. is the impurity concentration in the melt at the melting 

i, t 

interface, and v (p /p f ) is the mass velocity in the melt due to the 

C C I 

movement of the zone at rate v . Boundary condition (a) comes from 

c 

symmetry, (b) comes from the prohibition on transfer through the melt 
surface, and (c) and (d) represent material balances at the interfaces. 


B. Dimensionless Equations 

For convenience in solving the problem, equations were restated 
in the following dimensionless 'orms: 


where 


(1) Dimensionless Continuity Equation 


3v V 3v 
r , r . z 

W + T + Jz~ 


(A-13) 
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v a 
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and v “ y/p is the kinematic viscosity 
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(2) Dimensionless Momentum Equation 
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(A-15) 


The dimensionless boundary conditions for flow are: 


a) At Z - 2,/a » V - 0 , V - -V . 

r z cf 


b) At Z ■ -£/a , V « 0 , V - -V . or 

r z cf 


c) At Z - 0 , 


d) At R - 0 , 


») At R - 1 , 


3v 

V z " ® (only for V ^ ■ 0) 


3V 

V » ——2 

r 3 r 


3v 3fl 

v r - 0 , 55 - - **• s| 


(A-16) 


where 


'V p c 

r cf - \~tKz 


(A-17) 


is the dimensionless mass velocity in the melt due to movement of the 
zone, 

pa (T -T )a 
o m 


M' - 


(A-18) 
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is the dimensionless surface tension driving force parameter, and 


T-T 


8 - 


m 


T -T 
o m 


(A-19) 


The pressure terms in the momentum equations are eliminated 
by differentiating Equation (A-14) with respect to Z and Equation {A-15) 
with respect to R, and subtracting to obtain 


- ;p - ?) ■ * H ' (§ ' * $) * ‘•‘" 2 °> 


where 


3V 3V 


a) - s 


3Z 3R 


(A-21) 


is the dimensionless vorticity, 

<w‘ 


Gr 


h v 2 


is the Grashof number for temperature variations, 

ga(m )a 3 

Gr » ~2 


V 


is the Grashof number for composition variations, and 


(A-22) 


(A-23) 


m 

bl 


(A-24) 


is a dimensionless composition, a ■ (3p/3m)/p is the concentration 
densification coefficient, and m^ is the impurity concentration the feed 
solid would have if it had the same density as the melt. 




} 
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30 3c . 

The derivation of the terns, Gr^ ^ and Gr^ ^ in the 
momentum equation (A-20) is as follows: 

A 

The term - (a /V p) (3P/3R) in Equation (7~i4) is differentiated 
with respect to Z, and -(a^/V^) (g + ^ ) in Equation (A-15) i' dif- 

ferentiated with respect to R. Then, subtracting them by usi.>g 
Equations (a- 4) through (A-7) we obtain: 

3 


a 


V 


3p 
p 3Z 


+ 9 ) 



9p 

3z 


) 


. 2_ fl ?p'j + 8_ /I 
3z \p 3 r/ 3r \p 

= [_ 3 |l/ 3P o 3 p\I 

^2 3Z I p t 3R 3 r/| 

3 (1 3p\ 3 1 1 / ^ 3P \ 

3Z \p 3R / + 3R I p \ P f g 3Z/ 


V 


1 3 P 3lP 3_ /1\ 
p 3Z3R “ 3R 3z ^p ) ~ 


3 (T-T ) . -2 n 

0 m 1 3 P 

9P ^ T-% * 


3R 


p 3R3Z 


3P 3 /l\ 

3Z SR^pj 


gB (T-T)a' 
o m 


36. 

3R 


„ 36 

" Gr h 3 r 


8C 

The term Gr t— is generated by the same procedure by taking the 
in qR *" 

density change due to concentration gradient into account. 


t 


t - - 


The dimensionless momentum equation may alternatively be de- 
fined by a stream function such that continuity is automatically 
satisfied from the following relationships: 


,, 1 Ml 

V = — 7T- 
r R 9Z 


V = . iJJt 
z R 3R 

i / d 2 j> i 

R Vsb 2 R 


»_t 
3R 3 r 2 , 


From Eqs. (A-25) through (A-27) , the momentum equation (A-20) is 
expressed by a single dependent variable ^ as follows: 


i am / a 3 » i »5t + i!l. lift + 2_>t\ 

rSz U 2 3R r 3z 2 3r 3 r 3r 2 r 23 7 

1 ajp 3_ /a 2 » 3jf 1 1 3i|i 3 2 i|i 

- R 3r 32 ^ z 2 3r 2 ' R 3R / ' 3z ,*2 

1_ a / 3^, 1 ^ 

R 2 32 \ 3R 2 R 3 *l 

. L A . i JjL + + si _ l 

\ 3z 2 3r 2 r 32 2 3r 3z 4 3r 4 r 3r 3 


+ ? H ' ? Sr)* Gr h (3e/3Z) + Gr „, (3c/az) 

R 9R R / 


The boundary conditions are also transformed in terms of the 
stream function. Since the stream function along the axis is a 
constant, it may be taken as zero for simplification. On the basis 
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of this, the boundary c .ditions are 


b) At Z = -Jl/a » 

c) At Z = 0 , 

d) At R = 0 , 

e) At R = 1 , 


r 2 v 


a) At Z = Jl/a , 'I' = 


cf 


2 

R*V 


* = 


cf 


I 3^ n 
♦ *5z - 0 


♦-IS-o 


5z = 0 




or 


(for V _ = 0) 
ct 




_V 39 

2 ' Sr 2 3Z 


(A-29) 


(3) Dimensionless Heat Transfer Equations 

By introducing the previously defined dimensionless 
length and velocity parameters for r, z, v and v into Equation (a-9) 

Z JL 

we obtain 


30 

3Z 


+ V 


30_ 

3R 


i_ a 2 e , i 39 | a 2 e 

Pr 3 z 2 R aR 3 r 2 


(A-31) 


where Pr = C uA is the Prandtl number and 
P 

T-T 

m 


T -T 
o m 


or 


T-T 

E_ 

b T-T 
m c 


(A-19) 


(A-31) 
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The dimensionless boundary conditions for energy transfer ares 


a) At R = 0 , 


90 « 

9r = 0 


b) At R - 1 , |f - - (<8^T ni ) 4 -{0 c AT*T m ) 4 ] 


c) At 


i R “ 1 / 


Z * 0 , 


6 = • (line heat source at T ) 

AT o 


d) At Z = JL/a , 6 = 0 


e) At Z ■ -£/a , 9 = 0 , or 


f) At Z = 0 , = 0 (only for V ^ = 0) 


d*) At Z = b/a , 


* - rrtr- KQAT+T } 4 -{0 At+t } 4 
3z kAT I m c in j 

for radiation 


= 0 for insulation 


6 « 6 for infinite heat transfer 
coefficient 


where At is either (T -T ) if 0 is defined by Equation (A-19) , or 

o m 

(T ->T ) if 6 is defined by Equation (A-31) , and 9 [=(T -T ) /At] is 
me c c m 

the dimensionless temperature of the surroundings. 


( 4) Dimensionless Mass Transfer Equations 

The following dimensionless impurity concentration 


parameter C was employed with the previously defined length and 


i 



itfcaf 


velocity parameters for r, z, v , and v to transform Equacicn (A-ll) 

r z 


3c 

V ~ + V 
Z 3Z 


3c _ 1 _ 3^c 1 3c 

r 3R " Sc .2 R 3R 
3 r 


(A-33) 


where Sc - v/D is the Schmidt number. 

The dimensionless boundary conditions for mass transfer are 


a) At R - 0 , = 0 


b) At R = 1 , 3 r = 0 


c) At Z — -A/a , C 


i+l.b 

i,b = I P c\ I 

1 + ScV (k - — — JAz 
c\ ° P J 


d) At Z = Va , C. 


C ^Scv Az + c. . t 

etc l-l , t 


Scv Az(~ ) + l 
c \ p f/ 


where C ^ (= i s dimensionless impurity concentration in the 

feed crystal, subscript i denotes a Z grid point and V (= v a/V) is the 

c c 

dimensionless zone travel rate. The boundary conditions (A-35) (c) and 
(d) result from taking the finite difference equation form of the 
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boundary conditions (A-12) (c) and (d) , and expressing the impurity 

concentration in the crystal in terms of the equilibrium distribution 

coefficient k and the interfacial concentrations, 
o 

C. Computation 

One of the typical numerical procedures for calculating the 
stream function in an enclosed fluid is to solve for the vorticity 
function, and then find the stream function from vorticity by 
employing a successive iteration method (A- 5, A-6, A-7) . Rather than 
finding the stream function via the vorticity equation, a single 
step is necessary to compute the stream function when Equation (A- 28) 
is used. It is obvious that this procedure introduces equal or less 
computation errors them the former method since the errors caused 
in the process of computing the stream function via vorticity are 
entirely eliminated by using a single step computation. 

An approximate solution is obtained at a finite number of 
equidistant grid points having coordinates Z = (i-l)Az and 
R = (J-l)AR, where i and j are integers. The parabolic finite dif- 
ference equation was employed for Equation (A- 28) . The resultant 
computation molecule is shown in Figure A-2, in which twelve current 
values at the nearest neighboring grid pints are required for an 
explicit computation of the new stream function. Computations of the 
stream function at the grid points next to the liquid/solid boundaries 
and to the axis are possible due to the symmetry condition of the 
stream function at those boundaries. The stream functions at the 



m 



Figure A-2. The finite difference grid system and the 
computational molecules [ (I) for stream 
function, (II) for temperature and con- 
centration] used in the numerical calcula- 
tions. 



interior grid points next to the free liquid surface are calculated 
by the interpolation method , since the stream function at the free 
liquid surface is readily calculable using Equation (A-29) (d) . The 
Gauss-Seidel iteration method (A-12) was employed for computer calcu- 
lation. For M' < 350 we obtained convergent solutions when 11 x 11 
grid points were used. However, an increased number of grid pints 
were necessary for convergence with higher values of M', as shown in 
Table I. 

The dimensionless vorticities and velocities in the melt were 
computed by the finite difference equations derived from Equation 
(A-25) through (A-27) . The calculated flow velocities were utilized 
for computations of the temperature and impurity concentration fields. 

The same computational scheme was employed for the analysis of heat and 
mass transfer. The finite difference of Equation (A-30) for heat 
transfer, and of Equation (A-33) for mass transfer resulted in the com- 
putation molecule as shown in Figure A-2. The Newton-Raphson method 
(A-9) was used to iteratively solve for the temperatures at the surfaces 
because of the non-linearity of the thermal radiation boundary conditions. 
Equations (A-32) (b) and (d')> A rapidly converging solution for heat 
transfer was attained by using an overrelaxation parameter of 1.9 when 
the Prandtl number was low (metals) . In heat transfer analyses with 
high Prandtl numbers or in mass transfer analyses for high Schmidt 
number liquids (including silicon) , introduction of the overrelaxation 
parameter impaired the convergence had to be employed in order to 
obtain a convergent solution, at the expense of an increased number of 


iterations. 
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APPENDIX B 

ANALYSIS OF FREE CONVECTION 
ABOUT A SPHERE 


(Chong K. Kim) 

A. Dimensionless Equations 

The symbols were defined in Section III. Figure B-l shows the 
spherical coordinate system used for this problem along with its 

z 
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Figure B-2. Vector Analysis 
of Gravity on 
Sphere 


! 


orientation relative to the 
accelerational vector. 

Figure B-2 shows that the 
acceleration body force g 
may be decomposed into two 
components, g^ « -g cos 0 
and gg - g sin 0. 

The pressure in the 
momentum equation can be 
written as the sum of the 
pressure p' due to the fluid 
motion and that due to the 
fluid density far from the 
sphere p^, so that 
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and 


£e - 

3r 



. -g 


P» 


cos 6 







■ 9 P» 


sin 0 + 


&L 

r30 


(B-l) 


(B-2) 


The following dimensionless variables were employed: 

(2a) 3 ga(W o -Wj v ua .. va „ p'a 2 

Gr = - , Sc = — , u = — u = — p * 


>'a 2 


UV 


v a 

v -~ ■ p_ ? 

pV 


Finally, in terms of dimensionless variables the differential 
equations are: 

(a) Continuity ; 


9u 2 cot 0 1 3v 

9r r u v r 30 


= 0 


(B-3) 


(b) Momentum; 


3u „ 3U V3U V Gr . 

3? + u liTKSff ‘ 5 <T * cos 6 


-H 


3V 


3^ 4 3U 1 3^ . cot 0 3u . 2 „ 

17 717 ~ w 7 ° 

3V . V 3V UV Gr 


aT +O 3R + R30 + r"8* ,in0 

1 3P afv 2 3v 1 3^v 

"R30 3R 2 5 3R R 2 a0 2 

cot 8 3V 2 3U V 

R 2 39 R 2 00 ” R 2 sin 2 0 


(p -4) 


(B-5) 



134 


4 


(c) Mass Transfer; 



+ U 


W X V M 
9R + R 90 


1_ 

Sc 


9 2 <ft 2 . 9j>_ 8^4> cot 9 9$ 

9R 2 R 3r R 2 90 2 R 2 99 


<B-6) 


In order to solve the equations above we must know the appro- 
priate boundary conditions. Initially we have a completely stationary 
and homogeneous system. At the surface of the sphere W * W q . Along 
the Z-axis (9 = 0,ir), symmetry requires that the gradient of radial 
velocity and of concentration with respect to 0 cure zero. In mathe- 
matical form 


i) at t ■ 0 , 

4> = 0, U «* 0 , V - 0 

ii) at R ■ 1 , 

<j)«l, U-0, V » 0 

iii) at R = 10 

if-o , U-0, V-0 

iv) at 9 ■ 0 and it 
9d> 9u 

W " 0 » 3? " 0 ' V “° 


(B-7) 
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B. Method of Solution 

In order to eliminate the hydrodynamic pressure gradients, 

a k t* a 

r— and — r— • , from the momentum equations, they were differentiated with 
dR R du 

respect to 6 and R and subtracted from one another. 

We define the stream function and vorticity uj by 


U 


u 


1 

R 2 sin 8 


3<{» 

30 


i 


1 / 3U 3v 
R l 30 “ R 3 r 



-1 

R sin 8 3R 


(B-8) 


Since the changes in velocity and composition are expected to occur mainly 
near the sphere, we substitute Z ■ in R so that the pertinent equations 
become : 


Mass: 


24 + -z 

3t + e 

-2Z 

e 

Sc 


-2z 3 (e 2Z Ud>) 1 (Vd> sin 8) 
e 3Z sin 6 38 


+ 3| + lf| , cot 6 || 

3 z 2 32 36 2 30 


(B-9) 


Momentum : 


If- 




z 

z 

-Z 3 (Ue w) . 3 (Vu>) 




6 3 Z + 30 



3 Z 

3 2 (a)e Z ) 3(aie Z ) 

z 3 

1 3 (to sin 6) 


3 Z 2 _!r ' 

e 30 

sin 8 36 


-Gr 

8e Z 


cos 8 || + sin 8 || I 


(B-10) 


\ 
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Vorticity - Stream Function: 



(B— 11) 


Stream Function - Velocity: 


U 


1_ 1 
2Z sin t3 wJ 


V 


-1 1 3 ^ 

2Z sin 6 3Z 


(B-12) 


(B-13) 


Boundary Conditions: 

The vorticities at the solid-liquid interface (Z = 0) can 
be computed from a Taylor series expansion. Along the interface. 


iff [iSTT If ] and H are * ec °' “ ■ stnT ^2 • usi " 9 a Taylor 


expansion, 


9i|j, 


a 2 <l> 


2 a J <K 




i,m+l 




. Py l,m Az u *i,m (AZ) “ y i,m (AzT 
i,m 3Z I! az 2 21 3z 3 3! 

2 ^3.1 


S *i.n,2to . 3 \.» (2iZ1 . S3 *i.» (24Z) 3 , 

i,m+2 • *1,« + st tr * * —3 — + ••• 


At the surface Z - 0 , 4*. «* 0 and 


3\Jj. 

r x,m 

i,m 3z * 0 


We can derive the following equation 

3 \» »*l.^l-*1.^2 


(B-lfc) 


3Z 


2 (Az) 1 


In addition, the concentration and radial velocity at 9 ■ 0 and 


(B-14) 

(B-15) 


it are: 
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At Z » 0 , 


U) 


_ . n+1 , n+1 

n+1 _ ^i,m+l ^i,m+2 

i,m 2 (Az) 2 sin(iA0) 


(B-17) 


At 0 = 0 , 


.. „A n+1 A n+1 

.n +1 4$ .. .-9 . 

9 _ m+l,j m+ 2,1 

m f j “ 3 

4u n+1 .-U n+1 . 
n +1 _ m+l,j m +2 , j 

m,j 3 


(B-16) 


(B-19) 


At 0 = TT , 


. , n +1 .n +1 

4 a) . — © , 

, n +1 _ m-l,j m- 2 ,j 

*»* j ~ 3 


(B-20) 


4u n+ ^ ,-U n+1 


u n +l _ m-1, j m-2,j 


m,] 


(B-21) 


The computer program used was very time-consuming because we 
used a double-stability-checking DO loop. The convergence criterion 
used were as follows: 


MaX \ip n+1 -4> n . | 

Max 1 . n .1 


< -4 

= 10 


Max, n+1 n , 

< 10 " 3 

Max, n , 

j V i» j 


(B-22) 


(B-23) 



